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Non-axisymmetric oscillations of differentially rotating stars are studied using both slow rotation 
and Cowling approximation. The equilibrium stellar models are relativistic polytropes where dif- 
ferential rotation is described by the relativistic j-constant rotation law. The oscillation spectrum 
is studied versus three main parameters: the stellar compactness M / R, the degree of differential 
rotation A and the number of maximun couplings ^ max . It is shown that the rotational splitting 
of the non-axisymmetric modes is strongly enhached by increasing the compactness of the star and 
the degree of differential rotation. Finally, we investigate the relation between the fundamental 
^ . quadrupole mode and the corotation band of differentially rotating stars. 
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I. INTRODUCTION 



Differential rotation is believed to play an important role in nascent neutron stars [l[ as well as in binary neutron 
stars mergers 0, Until viscosity, turbulent motion and/or magnetic fields force the star to uniform rotation, 
dynamical or secular instabilities can be developed due to differential rotation. 
^Jy The study of the differential rotation phase of a neutron star life re-gained attention lately. Recent numerical simu- 
i— i , lations in Newtonian hydrodynamics 0> H B 0] show that in stars with high differentially rotating, non-axisymmetric 
dynamical instabilities can develop even for low values of f3 — T/\W\ ~ 0.01 — 0.08, where T is the rotational kinetic 
energy and W the gravitational binding energy. The so-called "low T/|TU|" instability can drive either one-armed 
^ I ■ spiral or bar-mode instabilities. The gravitational wave signal emitted via these instabilities might be detectable 
with the advanced generation of ground based detectors. If this instability develops in supermassive stars fp|, it 
Q\ , may produce gravitational waves detectable even by the interferometric space detector LISA. As suggested in Q , the 
low T/|TU| instability might be due to the shear instability that the corotating / mode develops when it enters the 
\^ , corotation band. A Newtonian study of this instability [iJJ, which is based on the analysis of the canonical angular 

■ momentum, confirms the presence of corotation modes, 
t^. | A key feature of the oscillation spectrum of a uniformly or differentially rotating star is the splitting of the eigen- 
modes, like the Zeeman effect in atomic physics. In Newtonian theory, the splitting of the non-axisymmetric oscillations 
of differentially rotating stars have been studied with perturbative techniques in [111 ] . A general relativistic description 
of the perturbations of differentially rotating neutron stars may alter, at least quantitatively, the Newtonian results. 
In general relativity, the perturbations of uniformly rotating stars are mainly studied in the so called slow rotation 
approximation. This approximation is practically valid for the study of all known pulsars even for those with rota- 
tional periods of 1.5-2ms. The slow rotation approximation fails to treat perturbations of neutron stars with periods 
below 1.5ms or fJ/f^Kcpicr > 0.25 — 0.3. Actually, there are no known pulsars with such small periods but it is not at 
all impossible that newly born neutron stars may have rotational periods below 1ms. 

The lack of perturbative studies for differentially rotating relativistic stars, the absence of any results for non- 
axisymmetric perturbations and the issue of low T/|TU| instability are the main motivations of this work. There only 
a few studies of fast rotating neutron stars using a perturbative approach [l2|, EH, Q EH] ■ While there is significant 
progress in the study of neutron star oscillations using evolutions of the non- linear equations [H, E3, 13 • Still the 
non-axisymmetric perturbations of differentially rotating relativistic stars have not been treated by any method yet, 
this paper is the first attempt to address the problem. Actually, in an earlier paper we derived, in the perturbative 
framework of general relativity, the equations describing the oscillations of a slowly and differentially rotating neutron 
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star in the Cowling approximation (from now on Paper I) . Here we study the effect of differential rotation on the 
oscillation spectrum. Moreover, by using the perturbative approach developed here we examine the relation between 
the low T/|W| instability with the existence of corotating modes. The main results of our study can be summarized 
in the following sentence: the rotational splitting of the non-axisymmetric modes is enhanced by stellar compactness 
and the degree of differential rotation. 

The structure of this paper is as follows. In section [II] we briefly describe the perturbative framework (more details 
are given in Paper I). In Section [Ml we address the boundary value problem and introduce the numerical techniques 
for solving it. In Section lTvl we present and discuss our numerical results derived for different stellar models. Section fVl 
is dedicated to the conclusions and to the possible extensions of this work. Finally, in the Appendix [XI we describe 
the structure of the eigenvalue equations, while the definition of the linear angular operators is given in Appendix [Bl 

As is common, throughout the paper we use geometrical units i.e. c = G = 1. Prime ( ' ) denotes derivatives with 
respect to the radial coordinate r and overdot ( ' ) denotes derivatives with respect to the time coordinate t. 



II. THE PERTURBATIVE FRAMEWORK 



Equilibrium configurations of differentially rotating relativistic stars have been already studied in the early '70s by 
Hartle [2(| and Will [2l[ . In the approach that we follow here, the background spacetime of a slowly and differentially 
rotating star assumed to be axially symmetric and can be described, at first order with respect to the angular velocity 
of the star, in Schwarzschild coordinates, by the following line element: 

ds 2 = -e 2v dt 2 + e 2X dr 2 -2ojr 2 sin 2 9dtd<f> + r 2 (d9 2 + sin 2 8d(f> 2 ) . (1) 

The scalar functions is, A depend only on the radial coordinate r and are determined by solving the Tolman- 
Oppenheimer-Volkoff (TOV) equations for a given equation of state (EoS). The metric function uj = uj (r, 9) describes 
the dragging of inertial frames due to stellar rotation and obeys the following ODE [2^ |: 



47r (e+p) re 2X - 



,.2A 



uj = -16ir(e + p) e 2X Q, 



(2) 



where e and p are the total energy density and the pressure respectively, I the harmonic index and SI = SI (r, 9) is the 
angular velocity of the star as measured by an observer at infinity. 

A differentially rotating stellar model can be constructed in the framework of the slow rotation approximation in 
two steps. First, one constructs the non-rotating stellar model by specifying the central energy density and the EoS 
and then solving the TOV equations. Afterwards, a law describing the differential rotation is specified, that is one 
assumes a specific functional form for Sl(r, 9) and then the metric function uj(r, 9) i s determined by solving ODE f2j). 
Here, we use the perturbative version of the relativistic j-constant rotation law (l9l l23j|: 



n(r,0) 



A 2 n 



v oj(r, 9)7 



A 2 



(3) 



where Sl c denotes the angular velocity on the rotation axis, while the parameter A specifies the degree of differential 
rotation of the star. For A — » co, the j-constant rotation law ([3]) leads to uniformly rotating configurations, i.e. 
SI — > Sl c . More details about the procedure of constructing equilibrium configurations are given in Paper I, where we 
have adopted a harmonic expansion of the variables uj and SI up to i = 3. 

The perturbation equations describing non-barotropic, non-axisymmetric oscillations of slowly and differentially 
rotating relativistic stars have already been derived in Paper I. We have actually used the so called Cowling approx- 
imation, i.e. the equations are derived by perturbing only the fluid variables in the energy momentum equations 
5 [T'' 1 *) = 0. The oscillations of the fluid are then described by five functions, namely the enthalpy H im , the three 



perturbed velocity components u\ 



(polar), 



(axial) and the radial component of Lagrangian displacement 



vector £ . The perturbation equations read: 
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where w = ft — ui . The explicit form of the linear angular operators are given in Appendix iBl 

As in the case of uniformly rotating stars Oil]!, Eqs. ((4][8|) form an infinitely coupled system of equations. In 
particular, differential rotation introduces extra couplings with respect to the uniformly rotating case. In the limit of 
uniform stellar rotation, f2 3 and o> 3 vanish and Eqs. (|4j{8j) reduce to the perturbative equations presented in [24[. In 
addition, like in the uniformly rotating case, these equations can be devided into two independent subsystems, the so 
called axial-led and polar-led [26|. The polar-led system is the one that includes polar perturbations with £ — \m\ + 2k, 
and axial perturbations with t = \m\ + 2k + 1, where k is an integer. Instead, the axial-led system is the one that 
includes axial perturbations with t = \m\ + 2k, and polar perturbations with £ = \m\ + 2k + 1. The overall parity 
of each system is preserved and is polar for the first and axial for the second. In the rest of the paper we will focus 
mainly on the polar-led perturbations and leave axial-led for another study. Finally, for barotropic oscillations, where 
the background adiabatic index T is equal to the adiabatic index of the perturbations Ti, the last equation is obsolete 
and the system reduces to four evolution equations. 

The study of the oscillations described by Eqs. ((4][8|) can be done either in the time domain as an initial value 
problem or in the frequency domain as an eigenvalue problem (Sec. IIIip . Our study was based on the last method 
although we have also tried numerical evolutions of the system. In this case, we implemented a numerical code based 
on the two step Lax Wendroff scheme [27j |. For some of the stellar models, the simulations show some numerical 
instabilities after an evolution time of ~ 20 — 30 ms. Before the developement of the instability, we were able to 
determine with a Fast Fourier Transformation (FFT) the mode frequencies which agree better than ~ 3 — 5% to the 
ones derived by an eigenvalue problem. 
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III. PERTURBATION EQUATIONS IN FREQUENCY DOMAIN 

The perturbation equations (|4j[8j) can be studied in the frequency domain by assuming a harmonic time dependence 
for the perturbation functions e~ lCTt , with a being the oscillation frequency. By replacing the time derivatives in 
Eqs. ([4J8]) by — ier (e.g. H — > —iaH) and transforming H tm — > iH e "\ Ug m — » iu| m and £ £m — > i^ m one obtains a 
purely real eigenvalue problem, which is formed by two ODEs for H lm and u{ m and three algebraic equations for u^ 1 , 
//.'" and 

In order to prescribe the eigenvalue problem in a more compact form we first define the following three vectors: 



(H tm ,u{ m ) T , s e,n = (ui m ,ui m ) T , v em =(£ em ,0) T , (9) 



and then three infinite dimensional vectors: 

u = (...y- im y m y +im ,...) T , s = (...,/- im ,/"V +i "\...) T , v = (...y- im ,y m y+ im ,.. f .(10) 

By using these definitions, the perturbative equations can be written in an operatorial form as follows: 

j} a ot -U + C-S + P-V, (11) 



dr 



S a -S = M-V, (12) 

. V = 7V-U, (13) 

where A^, C, T>, S a , M, TV, <2^ ot , are infinite dimensional linear operators which are defined in Appendix lAl In 
-4^ ot , S a , Q^ ot the subscript denotes the dependence of these matrices on the oscillation frequency a. These operators, 
as already mentioned above, couple perturbations with different harmonic indices £, leading to an infinite system of 
coupled equations. For a given azimuthal index m, the number of couplings in the system of equations (|1 lti!3|) can be 
controlled by the parameter £ max , which is the harmonic index £ where the tensor harmonic expansion of perturbations 
is truncated. In this way, the number of couplings is n c = f max — |m| + 1 as £ runs \m\ to £ m&K i.e. \m\ < i < £ max . 
The truncation of the couplings up to a certain £ max has been done for practical reasons since it is impossible to deal 
with infinity many couplings which do not contribute significantly in the final result. This approximation has been 
tested by studying the variation of the eigenfrequencies with respect to a varying value of ^ max - In general, we found 
that even by keeping only a small number of couplings the eigenfrequencies were converging (Sec. IIVB)l . 

One can form a boundary value problem by using Eqs. (|1HI13P together with an appropriate set of boundary 
conditions at the center and the surface of the star. In the parameter domain where the inverse of the operators S a 
and <2^ ot exist, both Eqs. (|12[) and (fT3)) can be solved for S and V respectively: 



S = 5- 1 -A^-U, V= (Q* *) 1 -TV-U, (14) 

and inserted into Eq. (jlip. Then the final matrix equation, together with the appropriate boundary conditions, can 
be used for the calculation of the eigenfrequency a. 

In the domain of frequencies where the two operators and Q^ ot are singular, Eqs. (fT2"|) and (fT"3|) cannot be solved 
for S and V. This issue will be discussed in more detail in Sec. IIII CI 



A. Boundary Conditions 

The perturbation equations (|1 ltil3|) can be solved as an eigenvalue problem for the variable a by fixing the boundary 
conditions at the center and the surface of the star. Regularity at the center of the star suggests that the various 
perturbation functions have the following behavior: 

H lm ~ r \ u[ m ~r l -\ u e 2 m ~r e , uf 1 - r l+1 . (15) 

For any combination of £ and m there is only one independent solution at the center, as H lm and are related via 
the relation: 

[£H im - {-a + mflt) u{ m ] | r=0 = . (16) 



5 



\m\ + 1 independent 



Moreover, since at the center of the star H lm and H l m are independent, there exist £ ma 
boundary conditions. 

The second boundary condition comes from the vanishing of the pressure's Lagrangian perturbation on the surface. 
This condition is satisfied when the perturbations obey the following ^ max — \m\ + 1 number of equations: 
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Notice that in the previous equation, the operator acting on the enthalpy perturbation H lm comes from the total 
derivative: 
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In fact, the last term of Eq. (fT8|) becomes the operator Cf 2 of Eq. fl7|) after the angular integration of the perturbation 
equations. 



B. Numerical Method 

The numerical technique that is used here has been already described earlier in [2~4| |. By selecting a certain value 
for £ max we allow n c = ^ max — |m| + 1 couplings between various harmonics where \m\ < £ < £ max . Then one must 
specify n c independent conditions at the origin for H lm that can be denoted by a vector 

H k = (H lm , H l ^ rn ) , where k = 1, n c , (19) 

and here we assume that Hi = (1, 0, . . . , 0), H2 — (0, 1, 0, . . . , 0) , . . . , TL nc — (0, . . . , 0, 1). Then the corresponding 
values for u\ m are estimated via Eq. (fP6|). By specifying a trial eigensolution cr, we integrate for each TCk Eqs. (|1HI13[) 
and we calculate the value of the Lagrangian pressure perturbation Ap e k on the surface. This quantity in general does 
not vanish and after n c integrations, we construct an n c x n c "pressure matrix" P. An eigenmode a is a solution of 
Eqs. (|1HI13[) that simultaneously satisfy Ap e = for any £. Since the n c initial values for Hk are independent, one 
can find appropriate coefficients such that 

fmax 

^2 a k &Pk = 1 for ever Y i = m... Cax ■ (20) 
k—m 

This is equivalent with the requirement that det P = 0. After integrating Eqs. (|1HI13P for a wide range of frequencies 
we isolate the zeros of det P = by a "root finding algorithm" . These zeros are the eigenmodes a of the problem. 

The eigenfunctions of the enthalpy and velocity perturbations can be constructed as a linear combinations of the n c 
solutions that correspond to the n c independent initial conditions Hk- The coefficients of these linear combinations 
are determined by solving the homogeneous system of linear equations described by Eq (|20[) . 



C. Continuous Spectrum 

For the range of frequencies a that the inverse operators of S a and Q^ ot do not exist, Eqs. (|12j|13|) are singular and 
consequently the eigenmodes cannot be determined using the frequency domain code. This is the same singularity 
that generates the continuous spectrum as in the r-mode studies [H, [23, |3(| HH, H2| , which have been carried out in 
slow rotation approximation. The continuous spectrum (CS) practically can be determined by solving the determinant 
of the 2n c x 2n c matrices that represent the operators S a and Q^ 3 ', for the mode frequency a. Let us first consider 
the operator S a for the case £ max = |m|. The determinant is a second order polynomial in a which has a double root 
of the form: 

2?Ti Til/ 

a = m (n x + 6O3) ~ £ , e + 1) fa + 6 u 3 ) - £ + (I, m, O3, ^3) , (21) 
where the function ^ is defined as follows: 

* EE -15 (3f> 3 - 4 W3 ) [tQ\ +X , m ~ {I + 1) Qlm] 

+ y n 3 [£ (£+!)-(£+ 1) {£ - 2) Q\ m - £ (£ + 3) Q 2 +1<m ] , (22) 
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Model 


p c MR 
(xl0~ 3 km- 2 ) (M Q ) (km) 


M/R 


(xlO^km" 1 ) 


BO 
CO 


0.662 1.400 14.151 
1.484 1.637 11.218 


0.146 
0.216 


2.700 
4.168 


TABLE I: The parameters of the background non-rotating stellar models of the B and C 
mass density, M the mass, R the radius and Q.k denotes the mass sheeding limit (Kepler 


sequence. Here p c 
frequency) . 


is the central rest 


flfc (xKT 3 ) 


ao ai <22 


as 




B 

C 


0.726 -1.272 3.972 
-2.239 1.811 5.388 


-0.853 
2.233 


3.970 
5.386 



TABLE II: The coefficients a,k of the Pade approximation of Eq (|27[) for both the B and C sequences of stellar modes. 



and the coefficient Qe m is given in Eq. (|B12[) . The interval of the CS is then cr|J < a < erf , where <x|J and erf are for 
the operator S a the values of a given by Eq. (|21~j) at stellar surface and centre respectively. In general, for £ max > m, 
there are £ max — |m| + 1 branches of the CS, which might as well overlap widening the band of it. Similarly to the 
uniformly rotating star [24| . for £ max — > oo the CS will cover all the spectrum. 

Similar behavior can be observed for the operator Q t ° t , where one can calculate a singular patch for ^ max = m via 
the following expression: 

a = m(n 1+ 6fl 3 ) - ymtt 3 [l - Q% m - Q 2 e+1 J . (23) 

Again in this case, the width of the CS is defined by the values of a given by Eq. ([23]) at the stellar center and surface 
i.e. cr^ < cr < crj? and the number of continuous spectrum patches when more coupling terms are included is again 
4iax - M + 1. 

It is obvious from the above analysis, that the operators S a and Q^ ot generate different bands of CS. For a 
barotropic EoS, Eq. (|13|) is trivial and the continuous spectrum is only generated from operator S a . Furthermore, 
when expression ([23"]) is satisfied the output of the angular integration of the total time derivative given by Eq. (fT8|) 
vanishes and the CS of the operator Q^ ot lies in the corotation band of oscillation modes. 

It should be noted that the existence of a continuous spectrum for uniformly rotating relativistic stars is highly 
debatable and it might be a result of the slow rotation approximation 32]. On the other hand, differential rotation 
is known to be associated with the presence of corotation points and existence of a continuous spectrum but this is a 
technical problem that may not be solvable within the limits of the slow rotation approximation. Thus we will not 
address the issues related to the presence of CS in the rest of this paper. 



Model 


T c 








J 


(A = 12 km) 


(ms) 


(xlO" 2 km" 1 ) 


(xlO -2 km" 1 ) 




(km 2 ) 


Bl 


1.719 


1.218 


0.435 


0.161 


0.935 


B3 


0.970 


2.160 


0.771 


0.286 


1.657 


B6 


0.657 


3.189 


1.139 


0.422 


2.447 


CI 


1.719 


1.218 


0.535 


0.129 


0.832 


C3 


0.970 


2.160 


0.949 


0.229 


1.474 


C6 


0.657 


3.189 


1.401 


0.338 


2.176 



TABLE III: Parameters describing rotating models of the B and C sequences, see also Table HI Here T c and Q c are respectively 
the period and the angular velocity on the rotational axis, while f2 e and e e represent the angular velocity and the dimensioless 
parameter, described in Eq. (|28)| . at the equator. The angular momentum of the star is denoted with J. 
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A[km] 

FIG. 1: This figure displays the dependence of <f> e = Q. e /Q. c on the parameter A for the B (solid line) and C (dashed line) 
models. 



Model BJ1 


T c 


n c 






J 


A (km) 


(ms) 


(xlO -2 km" 1 ) 


(xlO -2 km" 1 ) 




(km 2 ) 


5 


0.667 


3.143 


0.268 


0.099 


0.935 


15 


1.989 


1.054 


0.491 


0.182 


0.935 


50 


2.746 


0.763 


0.692 


0.256 


0.935 


100 


2.834 


0.739 


0.725 


0.268 


0.935 



TABLE IV: This table displays the main properties of the differentially rotating models BJ1 for four typical values of A. 



IV. NUMERICAL RESULTS 

In this section, we solve Eqs. (| 1 ltil3|) and investigate several aspects of the oscillation spectrum of non-axisymmetric 
oscillations of slowly and differentially rotating stars. In the first subsection we describe the procedure in constructing 
the background stellar models while in the second subsection the spectrum of non-axisymmetric oscillations is discussed 
and we provide tables with frequencies for several background configurations. 

A. Stellar Models 

In the slow rotation approximation the stellar models maintain their spherical shape and the rotation is treated 
as a small perturbation. This approximation is practically applicable to all known neutron stars even to those with 
rotation frequency of the a few hundred Hz. Here we specify a specific value for the central density for a given EoS 
and we generated a sequence of rotating models by varying the angular velocity of the star. For simplicity, we adopted 
the relativistic barotropic EoS: 

P = Kp r , e = p + ^ , (24) 

where p is the rest mass density and e the total energy density, while K and V are the polytropic parameters. We 
have chosen the so called B-models [HI , which represent differentially rotating polytropic stars with central rest 
mass density p c = 7.91 x 10 14 g cm~ 3 and polytropic parameters r = 2 and K = 217.86 km 2 . The nonrotating 
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m 


/ — 2 


^max — 3 


^max — 4 


^max — 5 




-2 


1.4793 


1.4766 


1.4768 


1.4769 


2 f 


-1 


1.7037 


1.7016 


1.7016 






1 


2.1271 


2.1252 


2.1246 






2 


2.2932 


2.2901 


2.2913 


2.2913 




-2 


3.5798 


3.5811 


3.5815 


3.5815 




-1 


3.8156 


3.8205 


3.8206 






1 


4.4221 


4.4278 


4.4273 






2 


4.6467 


4.6485 


4.6488 


4.6488 



TABLE V: Frequencies (a/2ir), in kHz, of the fundamental ( 2 f) and the first pressure mode ( 2 p 1 ) for the Bl stellar model with 
A = 12 km and for £ min < £ maK < ^ min + 3. 
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e 

FIG. 2: Frequencies (<r/27r), in kHz, of the 2 f and 2 p x nonradial modes with m = 1 for the B sequence of models having 
A = 12 km. The solid line corresponds to the purely first order slow rotation approximation, while the dashed line includes 
some second order corrections due to the component it 3 m of the perturbed fluid velocity. 




member of this sequence, which is denoted as B0, has the typical neutron star mass M = 1.4 Mq and radius 
R = 14.151 km. In order to investigate the dependence of the non-axisymmetric spectra on the compactness of the 
star, we constructed, in addition, a sequence of polytropic stellar models with 0.102 < M/R < 0.216. The more 
compact model (M/R — 0.216) of this sequence will be called C-model with p c = 2.0 x 10 15 g cm~ 3 . More details of 
the nonrotating members of the B and C sequences are provided in Table |TJ 

In a differentially rotating star, the angular velocity on the rotation axis Q c and the parameter A, which describes 
the degree of differential rotation, are the other two free parameters which need to be specified in constructing a 
sequence of rotating stellar model by using Eq. ((2|). The angular velocity at surface, is then related to the one at 
the axis by the following relation: 

n s = Q c <l>(A,6), (25) 

where 4> = <fi(A, 8) is a scalar function that depends on the law that describes the differential rotation. For the 
relativistic j-constant rotation law given by Eq. ([3]) this function reads: 



A 2 + e- 2 »r 2 sin 2 9 V ^ c 



4 +e rsiii S-^P . (26) 
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The form of 4> e on the equator is drawn in Fig. Q] for both B and C sequence of stellar models. These two curves, for 
practical reasons, can be Pade approximated by the following rational function: 

fit _ ao + a.A + ^A 2 
0e _ l + a 3 A + a 4 A2 ' [Z<) 

where the coefficients a,k arc listed in Table HT1 

The parameter that is commonly used in describing differentially rotating stars is the ratio T/|PF|, where T is the 
rotational kinetic energy and and W the gravitational potential energy respectively [331 ] . Another possible rotational 
parameter, which have been used in [34j|, is a function of the total angular momentum, the total mass and the two 
parameters defining polytrope. In the relativistic slow rotation approximation, the gravitational potential energy W 
can be accurately determined only at order O (f2 2 ) , where the monopole and quadrupole corrections of the gravitational 
mass and internal energy can be defined. Therefore, we define as dimensionless rotation parameter the following ratio: 

£e = nJ = ^ (A, 2 )=£c0( '2 } ' (28) 

where f2 e represents the angular rotation at the stellar equator and f2 k is the Kepler angular velocity that defines the 
mass shedding limit of a rotating star. In slowly rotating neutron stars (first order in f2), Qk can be approximately 
described by the angular velocity of a particle in a stable circular Keplerian orbit at the equator of a nonrotating star 
fix = \J M/R 3 . Finally, the quantity e c = Q c /Qk is the value of the dimensionless rotation parameter at the axis. 
The parameters of a few selected models, belonging to the B and C sequences, are shown in Table UTTl for A = 12 km. 
Differential rotation is described by a number of parameters i.e. tt c , J and A. The effect of these parameters on 
the oscillation spectrum will be studied for two families of background stellar models. In each family we fix one of 
the first two parameters and vary A. The first family of stellar models is the so called B-sequence in which we fix 
the angular velocity f2 c and we vary the parameter A. The properties of the various members of this family, e.g. 
their angular velocity Q e at the equator and of the dimensionless parameter e e , can be estimated from those of the 
A = 12 km model (Table [Hi]) by using the Pade expression ([27)) . The second sequence is constructed by keeping 
constant the angular momentum J in the B-sequence of models and by varying the parameter A. The models of this 
new family will be named BJ models. The angular momentum of the differentially rotating star can be determined 
by the following expression [201 ] : 

J = 2ir / / dr d0(e+p) e x ~ u w r 4 sin 3 6 . (29) 
Jo Jo 

After expanding the variable w = — lo in spherical harmonics up to index I — 3 and performing the angular 
integration, Eq. (|29j) becomes: 



8ir f R 

J=— j dr{e+p)e x -" vo x 
3 Jo 



(30) 



where only the i = 1 component of w contributes to the integral. In Table IIV1 the main parameters characterizing 
the BJl model are shown for selected values of A. 

Once the sequence of BJl model is known, the rotational velocity £lf ,Tn of another BJn model, with the same value 
of A, can be estimated via the following relation: 

^ c BJn - ^ C BJ1 ^ , (31) 
where J B1 and J Bn refer to the angular momentum of the Bl and Bn models respectively. 



B. Spectral Properties 

Any nonradial mode of a rotating star is characterized by its harmonic indices (£,m), where —£< m < £. The 
axisymmetric modes m = have been already discussed in Paper I, so here we focus on the non-axisymmetric 
oscillations (m / 0). 

The non-axisymmetric modes of a rotating star split in corotating and counterotating branches, whose pattern 
speed a v = ajm is respectively positive and negative. Note that the modes are assumed to behave as e ~ 1 ( <Tt - m </'). i n 
the slow rotation approximation, the eigenfrequencies of any (£, to) non-axisymmetric mode are linear functions of the 
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A = 12 km 




BO 


Bl 


B3 


B6 


Modes 


111 












-2 


1.883 


1.479 


1.171 


0.839 


2 f 


-1 


1.883 


1.672 


1.515 


1.349 




1 


1.883 


2.102 


2.276 


2.473 




2 


1.883 


2.293 


2.614 


2.970 




-2 


4.107 


3.579 


3.177 


2.744 




-1 


4.107 


3.808 


3.579 


3.334 




1 


4.107 


4.416 


4.657 


4.925 




2 


4.107 


4.647 


5.068 


5.534 



TABLE VI: Frequencies (u/2-k), in kHz, of the 2 f and 2 p x nonradial modes for selected B stellar models with A — 12 km. 
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FIG. 3: Frequencies (<j/27r), in kHz, of the 2 f and 2 p 1 nonradial modes for B stellar models with A = 12 km. 



rotation parameter e e or e c . Note that e c and e e are related through Eq. For a given stellar model the splitting 

of the modes, due to differential rotation, can be described by the following relation: 

a irn = a im ± a) e c = a^™ ± q(^, |m|, A) e e , (32) 

where a and a are scalar functions depending on the harmonic indices and the differential parameter A while do is 
the oscillation frequency of the nonrotating configuration with respect to an inertial observer. For rotational laws 
that are independent from the 9 coordinate (e.g. uniform rotation), relation (|32j) is simplified and becomes linear 
with respect to the azimuthal index m, 

a tm = ai m + ma i uni e c , (33) 

as it is known from the Newtonian theory (3o| . 

1. Dependence OTh Umax 

For any (I, to) nonradial mode, Eqs. (|1HI13[) can be integrated by fixing the parameters £ m i n and ^ max - The value of 
t-axin is fixed by the selection of index m (£ m i n = |to|), whereas the value of ^ max is chosen according to the maximum 
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m = 2 




I i I i I i I i I i I i I 

0,1 0,12 0,14 0.16 0,18 0,2 0,22 

M/R 

FIG. 4: The dependence of the splitting factor a on the stellar compactness is drawn both for the 2 f and the 2 p x modes. Each 
member of this sequence of polytropic stars has different compactness but the same fi e /£2 c ratio. 



number of desired couplings n c and it's upper limit is oo. Before proceeding in deriving any result it is important to 
test the dependence of the oscillation frequencies for a given value of £ on £ max . In Table fVl we present the results 
of such a study, we actually show how the values of the fundamental 2 f and pressure 2 pj modes of the Bl model 
depend on £ max . The variation of £ max has not affected the actual values of the modes by more than 0.2% — 0.3%. 
The same order of dependence of the mode frequencies on £ max has been also observed for the other models of the B 
and C-sequences. Therefore, £ — 2 mode frequencies can be accurately estimated by ignoring couplings with higher 
£'s i.e. by keeping £ max = 2. 

In paper I, we discussed the role of the coupling between the polar and axial perturbation functions in the axisym- 
metric case. Actually, the component of the velocity perturbation Wg™ is intrinsically of O (SI) order, (see Eq. 0). 
This implies that there are "hidden" second order in f2 terms which induce some extra couplings into Eqs. (|5l6p 
that can influence the results. The same type of problem is also present for the non-axisymmetric perturbation, 
and will especially influence specific modes as for example the (£,m) — (2,1). Thus we studied the (£,m) = (2,1) 
eigenfrequencies either by using the full coupled system of equations, i.e. (£ m i n , £ max ) = (1,2), or by removing the 
"hidden" second order couplings, i.e. (f m in,^max) = (2,2). The results are shown in Fig. and the similarities to the 
axisymmetric results of Paper I are obvious. That is, as we expected the dependence of the mode frequencies on the 
rotation rate is not any more linear due to the presence of these implicit second order rotational terms. This effect is 
actually more pronounced for the 2 f mode. 

We have chosen to neglect all of these "hidden" second order terms in order to get the expected linear relation 
between the oscillation and the rotational frequencies. 

2. Eigenfrequencies 

Here, we study how the mode frequencies depend on the rotation parameter e e , the compactness M/R of the star 
and the degree of differential rotation A. 

Let us start by choosing the B sequence of stellar models for a fixed value of A = 12 km, these are the equilibrium 
configurations already used in [TtI [l8j . Due to rotation, these modes are split in two symmetrical branches charac- 
terized by the azimuthal number m. In Fig. [3] the splitting induced by rotation for the £ = 2 fundamental and first 
pressure modes is shown while the actual values of the eigenfrequencies for these sequence of models are given in 
Table ED 

Stellar compactness affects significantly the mode splitting, but one should be careful on how to quantify this effect. 
In fact, the angular velocity profile of the stellar model depends on the compactness and it is inconsistent to create a 
sequence of models by just keeping A constant. Instead, we have choosen to compare stellar models according to the 
ratio between the angular velocity at the equator and at the rotational axis, i.e. by keeping <p e = Sl e /fl c constant. 
We then consider the B sequence of stellar models with A = 12 km as reference for generating other polytropic stars 
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described by the Eq. (|24| . All these models have <j) e = 0.3572, which for the C models corresponds to A = 10.08 km. 
In Figure S] we show the dependence of the "splitting" coefficient a (|3"2"|) on the stellar compactness for the 2 f and 2 p : 
modes. It is obvious that the splitting is enhanced by the compactness and it can be easily be twice as large for very 
compact neutron stars. In principle, within the gravitational wave asteroseismology [36], this effect of rotation can 
be used, to infer the rotation and compactness of the oscillating neutron star. 

Finally, we investigate the dependence of the non-axisymmetric modes on the parameter A describing the degree 
of differential rotation. We consider the two sequences of stellar models described earlier in Sec. IIV Al These are, 
the B sequence of stellar models with tt c constant and the BJ sequence with J constant. In Table IVII1 we report 
the frequencies of the 2 f and 2 p 1 modes for some of these models. In Fig. [5l it is shown as the value of A affects 
the splitting of the 2 f and 2 p 1 modes. The splitting of the eigenfrequencies for stellar models with high degree of 
differential rotation depends strongly on the rotation parameter e e . Furthermore, it is worth noticing that the B and 
BJ models show the same dependence of a with respect to A. This behavior is expected since the differential rotation 
is described in both stellar sequences by the relativistic j-constant rotation law ([3]) and the no-rotating model of the 
sequences is the same. 

As it has been mentioned in the introduction, differential rotation plays an important role on the onset of rotational 
instabilities. Non-axisymmetric pulsations of rotating stars can become dynamical or secular unstable and enhance 
gravitational wave emission, e.g. see [371 [H| and references there in. These instabilities appear when the rotational 
parameter fj — T/\W\, reaches some critical value (3 C . The onset of secular and dynamical bar mode instabilities 
in Newtonian incompressible and uniformly rotating bodies is at j3 c = 0.14 and j3 c — 0.27 respectively. The secular 
instabilities are driven by dissipative processes, such as gravitational radiation via the Chandrasekhar-Fricdman-Schutz 
(CFS) mechanism and viscosity. 

Actually, the so called "low T/W" dynamical instability appears in stellar models with high degree of differential 
rotation. Several studies carried out with nonlinear hydrodynamical codes [U, 0] and perturbative methods 0, [13] 
suggest a strong correlation between the onset of the low T/W instability and the presence of corotation modes. By 
definition, these modes have their pattern speed equal to the local angular velocity of the star, i.e. a/ra — fl(r,9). 
For a differentially rotating star the possible corotation band of the spectrum is given by the interval Q e < a/m < il c , 
which is obviously larger for highly differentially rotating stars, i.e. with smaller A. For the B sequence of stellar 
models, wc estimated the required differential and rotation parameters for having corotation modes. Initially, we 
show in Fig. [5] the cases A — 5 km and A — 12 km. As expected, for A = 5 km the £ = m = 2 fundamental mode 
"corotates" for smaller rotation rate e e than the A = 12 km configuration. Fig.[7]shows the variation of the £ = m = 2 
fundamental mode with respect to the parameter A. All these models with the exception of the Bl and B2 present a 
corotation 2 f mode. The upper limit of the parameter A for having corotation is indicated with A c and it is illustrated 
with a circle in Fig. [7] Therefore, the star has corotation modes when A < A c . In Table IVIIII the critical values 
A c and the corresponding eigenfrequencies are listed for various models of the B sequence. When non-axisymmetric 
modes are into the corotation band the eigenvalue problem is mathematically singular. Still, for some stellar models 
it was possible to isolate the eigenmode frequencies by carefully studying their eigenfunctions into the continuous 
spectrum |10| . 

V. CONCLUSIONS AND DISCUSSION 

We presented a first comprehensive study of non-axisymmetric oscillations of slowly and differentially rotating 
neutron stars in the perturbative framework of General Relativity. By using Cowling approximation, we examined the 
spectral properties of polytropic stars and investigated their dependence on four main parameters: stellar compactness 
M / R, stellar rotation rate at the equator e e , degree of differential rotation A and the maximum number of perturbative 
couplings ^ max . In accordance with the first order slow rotation approximation, the non-axisymmetric modes exhibit 
a linear splitting with respect to the rotational parameter e e . However, in some of the eigenmode patterns appears 
a quadratic deviation with respect to the expected linear behavior. This is due to the presence of "implicit" second 
order perturbative terms in the perturbation equations. We have identified and then neglected these terms in order 
to be consistent with the order of approximation adopted for the background spacetime. Moreover, we show that the 
non-axisymmetric spectrum can be described by including only a small number of couplings between perturbation 
functions. For instance, we determined the quadrupolar spectrum with an accuracy to better than 1% by setting 
f max = 2, which corresponds to the lowest possible number of coupling terms in the perturbation equations. 

We found that both differential rotation and stellar compactness affect the non-axisymmetric spectrum. In fact, the 
rotational splitting of the non-axisymmetric modes is enhanced by stellar compactness and the degree of differential 
rotation. 

For the study of the "low T/W" instability we calculated the corotation band of some polytropic models and the 
necessary rotational configuration for having a corotating quadrupolar fundamental mode. Moreover, we verified 
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Modes 


m 


A (km) 


Bl 


B3 


B6 


BJ1 


BJ3 


BJ6 




1 


5 


1.965 


2.029 


2.101 


2.097 


2.269 




2 f 


1 


12 


2.102 


2.276 


2.473 


2.102 


2.276 


2.473 




1 


50 


2.244 


2.535 


2.864* 


2.107 


2.285 


2.484 




1 


100 


2.256 


2.101 


2.896* 


2.107 


2.285 


2.483 




2 


5 


2.027 


2.138 


2.627 


2.274 






2 f 


2 


12 


2.293 


2.614 


2.970 


2.293 


2.614 






2 


50 


2.615 


3.220 


3.917* 


2.334 


2.698 


3.113 




2 


100 


2.647 


3.283 


4.018* 


2.338 


2.707 


3.128 




1 


5 


4.226 


4.317 


4.416 


4.412 


4.648 


4.908 




1 


12 


4.416 


4.657 


4.925 


4.416 


4.657 


4.925 




1 


50 


4.609 


5.011 


5.465* 


4.419 


4.665 


4.940 




1 


100 


4.624 


5.038 


5.508* 


4.418 


4.664 


4.938 




2 


5 


4.295 


4.438 


4.593 


4.586 


4.949 


5.486 




2 


12 


4.647 


5.068 


5.534 


4.647 


5.068 


5.534 




2 


50 


5.114 


5.947 


6.912* 


4.728 


5.229 


5.799 




2 


100 


5.164 


6.044 


7.067* 


4.737 


5.247 


5.829 



TABLE VII: The mode frequencies (a/2n) of the fundamental ( 2 f) and the first pressure mode ( 2 p x ) in kHz, for the B and BJ 
sequences of stellar models. The B6 models for A > 38 km are rotating faster than the mass sheeding limit, their eigenmodes 
are labeled with a star. In the BJ3 and BJ6 models, some of the eigenfrequencies cannot be determined as they are inside the 
continuous spectrum. Therefore, we leave a blank space in the table. 




A [km] 

FIG. 5: The dependence of the splitting factor a on the differential rotation parameter A for the B and BJ models is drawn for 
the 2 f and 2 p l modes. The solid and the dashes lines represents respectively the 2 / and 2 pi modes of the B models, whereas 
the open circles and triangles the 2 / and 2 pi modes of the BJ models. 



the results of Newtonian studies, i.e. that the value of the rotational parameter e e , which is required for having a 
corotating fundamental mode, is inversely proportional to the degree of differential rotation. The onset and the details 
of the low r/IM^I dynamical instability of differentially rotating stars will be the subject of a future investigation. 
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FIG. 6: The upper panel displays the i — m = 2 f-mode (solid line) for the B sequence of stellar models with a high degree of 
differential rotation A = 5 km. The limiting lines of the corotation band are represented in dashed lines. The lower panel is 
similar but for A = 12 km. 



Models 




A c 


a/2% 






(km) 


(kHz) 


B3 


0.042 


3.898 


2.0615 


B4 


0.162 


7.422 


2.4105 


B5 


0.302 


10.182 


2.7336 


B6 


0.446 


12.534 


3.0439 


B7 


0.610 


15.095 


3.3454 


B8 


0.778 


17.625 


3.6707 



TABLE VIII: Critical values of the differential rotational parameters for having corotation modes. These values correspond to 
the I — m — 2 f-mode. 
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APPENDIX A: COEFFICIENT MATRICES 



In Sec. IIII1 we wrote the fluid perturbations as three infinite dimensional vectors U, S and V, see Eq. (fT0|) . These 
vectors can be written in terms of the 2- vectors u lm , s lm and v lm that are defined in Eq. ([9]). They contain the fluid 
perturbations with harmonic indices (I, m) . These definitions enable us to describe the perturbation equations in a 
more compact form by introducing a set of infinite dimensional linear matrix operators, i.e. A 1 " 1 , C, T>, S a , A4,Af, Q^ ot . 
In this Appendix, we provide the transformation laws of these operators. Actually, it is much more convenient to 
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FIG. 7: For B stellar models, the figure displays the variation of the I — m — 2 fundamental mode with respect to the 
differential parameter A. In any sequence, the circles denotes when the fundamental mode goes in corotation. For A lower 
then these values, the fundamental mode is always in corotation. The dashed line denotes instead the stellar models that are 
rotating faster than the mass shedding limit. 



define how these operators act on the single vectors i.e. the ones defined in Eq. ©, rather than dealing with the 
infinite dimensional vectors U, S and V given by Eq. (| 10[) . Let us start with the operator _4^ ot used in Eq. (fTTj) , which 
can be written as the sum of two operators A a and B ±2 defined as follows: 



A a = AI, 



B 



±2 _ 



±2 



where A is the following 2x2 matrix: 
/ 



1)1/ c 



'„-2 



V 



(a — m (Qi + -3 + 2m (w\ + 6073) 



e a(A-iO _2 +A '-2z/' + 4 



(Al) 



(A2) 



where I is the identity matrix and C] 1 is an operator defined in Eq. (|B5[) . The action of A a on u im is given by 



A a u lm 1 — ► Au im . 



whereas by the definition of Cf , the operator can transform u in three ways: 



B ±2 u em 1— » B [-Qe-! 



e-2m 



> (1 Qlm Qe+lm) 



where the coefficient Qe m is defined in Eq. (|B12j) and B is the following 2x2 matrix: 



B 







(A3) 



(A4) 



(A5) 



Note that the operator B ±2 couple perturbations with different index I while A a does not. 
The other two operators C and V of Eq. (fTTj) are defined as follows: 



C = Col © Ci/lf 1 © c 2 cf 2 



±3 



D J, 



(A6) 
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where Cf 1 , Cf 2 and Cf 3 are given in Sec. [Bj and the matrices Co, Ci, C 2 , C 3 , Dp have the following expressions: 



C 



-m(/x + 6/ 3 ) 
£(£+l)e 2X r- 2 



f mh 




c, = 



Ci 



f h 




-(/1 + 6/3) 




D 



v'c- 2 




and 



h = 2 - - v' )w x - wj , 



/ 3 = 2 - - v' w 3 - w 3 



The operation of for j = 1, 2, 3 on u tm follows the previous description as for the operator B 



±2 



In Eq. (|12|) the operators S a and Ai are defined as follows: 



S a = E CT Z © Ki£ 



±1 



JiCf 2 



J 2 £ 3 b2 



±3 
2 ; 



X = M J © Mi^ 1 © M 2 £± 2 © M 3 £ ±3 



3 ' 



where the quantities with a hat are matrices. Those associated with the operator S a have the following form: 



E„ = 



2 771 

-<7 + TO (ill + 6O3) — (lU! + 6a7 3 ) 



l 2x2 , 



2 (gTl + 6n7 3 ) / 1 
4 (4 + 1) 1 1 



K 2 = 



15to 2 / 2tn 3 + n 3 
1(1+1) [ 2w 3 



(A7) 
(A8) 

(A9) 

(A10) 
(All) 

(A12) 

(A13) 



Ji 



15to / Q, 3 - 2cj 3 
i{£ + l) [ 2w 3 



15to / 2w 3 
l{t+l)\ Q 3 - 2uj 3 



(A14) 



'3 = 



15 m 
"T^(£ + 1) 



n 3 i 



3A2X2 , 



w = 



15 



n 3 - 2u 3 



£(£ + 1) \ 2w 3 



here I 2x2 is the identity matrix of rank 2. While the matrices associated with the operator Ai are: 



Mo 



1 



t(l+l) mr 2 e- 2A (g 1 +6g 3 ) 



1(1 + 1) \ 







Mi = 



1 











£{£ + !) \ r 2 e- 2X ( 9l + 6g 3 ) 



M 2 



to / -fr 2 e- 2X g 3 



£(l+l)\0 



M, 







£(£+1) \ 



fr 2 e" 2 \ 3 



where 



gi = 2 ( 1/ ) tui + w\ , g 3 = 2 ( v' \w 3 + w' z . 



J" J \r 
Finally, we consider the operators Q^ ot and Af of Eq. (|T3J), which exists only in the non-barotropic case: 

= Q a I © Qi£f 2 , 
N = N X, 

where the matrices Q a , Qi and Nq have the form: 



Qn 



-a + m(fli +6n 3 ) 




Qi 



^TO^ 3 













e 2 ^ 2A (1 - ^) v' 



(A15) 

(A16) 
(A17) 

(A18) 



(A19) 
(A20) 



(A21) 



However, the definition of the vector v in ^ and the values of the matrix coefficients in (|A21[) lead to simpler 
equations: 
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-a + to (fli + 6O3) — — to n 3 cf 



±2 



(A22) 
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APPENDIX B: ANGULAR OPERATORS 

In this Appendix, we present the definition of the linear an gul ar operators with i,j £ N, which are used in the 
angular integration of the perturbation equations derived in [19( and couple the various perturbation functions. For 
convenience we don't write down the detailed expression of each operator, but rather the final part that shows the 
couplings that each one of them introduces on a generic harmonic component perturbation Ai m . 

The operators that introduce couplings of the form t ± 1 are: 

£\ X Ai m = (£—\) Qim^e-im — {£ + 2) Qe+lm^e+lm j (Bl) 

C 2 X Al m = — {£ + 1) Qe m Ai-l m + IQi+lmAl+im , (B2) 

C^Aim = (£-l)(e + l)Q im Ai- lm + £(£ + 2)Q i+lm A i+lm , (B3) 

C^Aim = QtrnAe-lm + Qe+lmAe+lm , (B4) 

while the following operators introduce couplings of the form t ± 2, 

£\' ''Aim — —Qe-lmQimAi—2m + (l — Qim ~ Qt+lm) -Atm ~ Q ' imQ 'f+lm^+2m i (B5) 

£-f 2 Al m = — {£ + 1) Qi-lm,Qim-Ae-2m + [f-Qj+lm — + 1) Qim] -^lm + £Q t+lmQ t+2mAi+2m , (B6) 

£3 2 Al m = (£ — 2) Ql-\ m QlmAt-2m + [^Qi+lm — + 1) Qim] Aim ~ [£ + 3) Qe +lmQ I +2mAi+2m , (B7) 

£± 2 Al m = — {I — 2) {I + 1) Ql-lmQtrnAe-2m — P- {£ + 3) Qi + i m Qe+2mAe+2m 

+ [£ (£ + 1) - {£ + 1) {£ - 2) QL - £ (£ + 3) Q, 2 +lm ] ^ m . (B8) 
Finally, the following operators introduce couplings of the form £ ± 3, 

— 2mQi— ImQirn Ag~Zm + + 4) Q g+lmQ £+2mQ i+3m A(.-\-% m 
+Qlm [tQi-lm + - 1) l 1 _ Qin - Q&flm)] A-lm 

-Qmm [(* + 1) Q, 2 +2m + (* + 2) (l - QL - Qf +lm )] Al + x m , (B9) 

^^(m = — (£ — 3) (£ + 1) Qe-2mQe-lmQimAi-3m ~ £ {£ + 4) Qi+lm Q^+2mQ^+3mA+3m 
+ Qftn [* (* + 1) QLlm " (* " 1) (* + 1) OL + * (t - 1) Q'+lm] A-lm 

+Qe+im [(£ + 1)(£ + 2) QL - 1 (i + 2) QI+ Xm + * (i + 1) Q, 2 +2m ] A+im , (BIO) 

Cf 3 Ai m = {£ + 1) Qi— 2mQ^-lmQ^mA-3m — ^Qi+lmQi+2mQi+3m A+3m 
-Q^m [(* + 1) + ^ +lm - (* + 1) {Qi-lm + Qim)] A-lm 

+Q^+i m [* + + 1) QL - £ (Q| +lm + Q| +2TO )] A+im , (Bll) 

where Qf m is defined as 



_ (£-m) (l + m) 

V (2€-l)(2€+l) " (B12) 



A useful relation between operators is the following 
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